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ABSTRACT 

Context. Internal gravity waves (hereafter IGW) are known as one of the candidates for explaining the angular velocity profile in the 
Sun and in solar-type main-sequence and evolved stars, due to their role in the transport of angular momentum. Our bringing concerns 
critical layers, a process poorly explored in stellar physics, defined as the location where the local relative frequency of a given wave 
to the rotational frequency of the fluid tends to zero (/. e. that corresponds to co-rotation resonances). 

Aims. IGW propagate through stably- stratified radiative regions, where they extract or deposit angular momentum through two pro- 
cesses: radiative and viscous dampings and critical layers. Our goal is to obtain a complete picture of the eff'ects of this latters. 
Mettiods. First, we expose a mathematical resolution of the equation of propagation for IGWs in adiabatic and non-adiabatic cases 
near critical layers. Then, the use of a dynamical stellar evolution code, which treats the secular transport of angular momentum, 
allows us to apply these results to the case of a solar-like star. 

Results. The analysis reveals two cases depending on the value of the Richardson number at critical layers: a stable one, where IGWs 
are attenuated as they pass through a critical level, and an unstable turbulent case where they can be reflected/transmitted by the crit- 
ical level with a coefficient larger than one. Such over-reflection/transmission can have strong implications on our vision of angular 
momentum transport in stellar interiors. 

Conclusions. This paper highlights the existence of two regimes defining the interaction between an IGW and a critical layer. An 
application exposes the effect of the first regime, showing a strengthening of the damping of the wave. Moreover, this work opens new 
ways concerning the coupling between IGWs and shear instabilities in stellar interiors. 
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1. Introduction 

Thanks to helio- and asteroseismology, we are able to extract a 
huge amount of informations about solar and stellar structures 
and compositions (e.g. Turck-Chieze & Couvidat (2011); Aerts 
et al. (2010)) and their internal diflPerential rotation profile 
(Garcia et al. 2007; Beck et al. 2012; Deheuvels et al. 2012). We 
know that internal rotation modifies the stellar structure since 
it generates flows, instabilities and chemical elements mixing, 
which modify stars evolution, for example their lifetime and 
their nucleosynthetic properties (e.g. Maeder 2009, and refer- 
ences therein). Moreover, in order to understand the obtained 
solar/stellar rotation profiles and the related rotational history, 
it is essential to develop a complete theory incorporating the 
diflTerent angular momentum transport mechanisms occuring in 
stellar interiors (e.g. Mathis 2010, and references therein). In 
this context, one legitimely asks the question about the origin 
of the internal rotation profiles. Four main processes are re- 
sponsible, in diflferent ways, for the secular angular momentum 
transport in radiative interiors. First, a large-scale meridional 
circulation is driven by structural adjustments of stars, external 
applied torques and internal stresses (e.g. Zahn 1992; Mathis 
& Zahn 2004; Decressin et al. 2009). Next, rotation profiles 
may be subject to diflTerent hydrodynamical shear and baroclinic 
instabilities and turbulence (e.g. Knobloch & Spruit 1982; 
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Talon & Zahn 1997; Maeder 2003; Mathis et al. 2004). Then, 
fossil magnetic fields, trapped during early-phases of stellar 
evolution once radiation zones have been formed (Braithwaite 
& Spruit 2004; Duez & Mathis 2010), can transport angular 
momentum through large-scale torques and Maxwell stresses 
(e.g. Gough & Mclntyre 1998; Mathis & Zahn 2005; uaraud 
& Garaud 2008; Strugarek et al. 2011). Finally, IGWs excited 
at the convection/radiation boundaries constitute the fourth 
mechanism able to transport angular momentum over large 
distances in stellar radiation zones (e.g. Press 1981; Goldreich 
& Nicholson 1989; Schatzman 1993; Zahn et al. 1997; Talon 
& Charbonnel 2005; Mathis & de Brye 2012). Note that all 
these processes are not necessarily present at the same time 
everywhere in the H.-R. diagram. Moreover, they act on various 
characteristic timescales in stars of diflTerent masses and ages. 

The object of this paper is the propagation of IGWs and the 
way they interact with the shear (i.e. the diflferential rotation) 
of the surrounding fluid. They are common in the terrestrial 
atmosphere and oceans (Eckart 1961; Chapman & Lindzen 
1970), that is why they are pretty well known in Geophysics. 
We will use this advantage for their study in the stellar case. We 
here draw attention to the mechanism whereby IGWs exchange 
energy with the mean flow, independently from other dissipative 
processes such as thermal and viscous diflTusion. Indeed, when 
the frequency of excited waves is of the same order as the 
angular velocity of the fluid (we will see the accurate definition 
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later), a phenomenon of resonance occurs, which affects the 
properties of both the wave and the shear of the surrounding 
fluid. This phenomenon is called a critical layer. Under the 
assumption of a perfect fluid (neither heat conductor nor 
viscous), Booker & Bretherton (19o/) and Lindzen & Barker 
(1985) have provided first results about critical layers in the 
geophysical case. They have shown that, depending on the 
value of the Richardson number of the fluid, which compares 
the relative strength of the shear and of the stable stratification 
(see Eq. 3), the waves might be either attenuated or reflected by 
the critical layer. This reflection may even be an over-reflection 
together with an over-transmission (see also Sutherland & 
Yewchuk (2004) for a laboratory evidence). In the same time, 
Koppel (1964); Hazel (1967); Baldwin & Roberts (1970); Van 
Duin & Kelder (1986) have completed this work in taking into 
account the conduction of heat and the viscosity of the fluid. 
Surprinsingly, their conclusions about the role of critical layers 
are identical. However, all these authors have produced their 
study in cartesian coordinates, assuming that the stifl'ness of 
the domain was thin enough to neglect its curvature. In the 
case of stellar radiation zones where critical layers may play 
an important role (e.g. Barker & Ogilvie 2010; Barker 2011; 
Rogers et al. 2012), these equations should thus be generalized 
to the case of spherical coordinates to be able to treat deep 
spherical shells. 

Therefore, after exposing our notations and assumptions (§2.), 
we present a complete mathematical study of critical layers in 
the case of a perfect fluid (§3.) and of a non-perfect fluid (§4.) in 
spherical geometry. Then, we determine the related mean verti- 
cal flux of angular momentum transported by IGWs (§5.). Next, 
we implement our theoretical results in the dynamical stellar 
evolution code STAREVOL and we apply our formalism to the 
evolution of a one solar-mass star (§6.). Finally, we present the 
conclusion and perspectives of this work (§7.). 



2. Definition, notations and hypotheses 

The star we consider is composed (at least) with a convective 
and a radiative region. IGWs are excited at the boundary 
between these two regions. They propagate in the radiative zone 
and are evanescent in the convective zone (e.g. Press 1 Q'^ ^ ). As 
explained in the introduction, in Geophysics, the study of IGWs 
is usually made in cartesian coordinates, justified by a local 
approach. In the stellar case however, a global approach using 
the spherical coordinates (r,6,(fi) is necessary. 

In the frame of Zahn (1992), we choose a shelular angular 
velocity for the studied star's radiation zone: ^(r,0) = ^(r), 
considering that, because of the strong stable stratification, the 
shear instability decreases the horizontal gradient of the angular 
velocity (e.g. Talon & Zahn 1997; Maeder 2003; Mathis et al. 
200 ), which consequently can be considered as only dependent 
of radius. For the moment, we neglect the action of Coriolis 
and centrifugal accelerations while the Doppler shift due to 
diff'erential rotation is retained. We also neglect the action of a 
potential magnetic field. 

Now, we need to introduce some quantities to describe the prop- 
erties of the fluid in the studied radiative zone. Each scalar field 
X is written as 



where we have introduced its horizontal average on an isobar, X, 
and its associated fluctuation X'. The thermodynamic variables 
employed are the density p, the pressure p, the temperature T 
and the specific entropy S . Next, the stratification is described in 
terms of the Brunt- Vaisala frequency 



p dr 



1 dp 
Yip dr 



(2) 



where g{r) is the mean gravity in the Cowling approximation 
where fluctuations of the gravific potential are neglected (see 



Cowling 1941), and Yi 



d\np\ 

— — r the adiabatic exponent. The 
amp ' 



relative importance of the stable stratification restoring force 
and the shear destabilizing eff'ects is quantified thanks to the 
Richardson number 

Ri=^^. (3) 



X (r, 6>, ipj) = X (r) -h X' (r, 6, ip, t) , 



(1) 



dr 

When Ri is small, the velocity shear overcomes the stabilizing 
buoyancy and turbulence and mixing occur (e.g. ^alon & Zahn 
). On the contrary, when Ri is large, the fluid remains 
stable. Finally, for the case of a non-perfect fluid (§4.), we 
introduce the viscosity v of the fluid and the coefficient of 
thermal conductivity k. These notations will be recall at the 
proper moment. 

IGWs themselves are characterized with their relative frequency 

cr(r) = (Tw + mAQ(r), (4) 

where (Tw is their excitation frequency (from the base of the 
convective zone in low-mass stars or the top of the convective 
core in intermediate and high mass stars), m corresponds to a 
Fourier expansion along the longitudinal direction (c.f. Eq 10) 
and AQ(r) = Q(r) - Qcz is the difference between the angular 
velocity at the level r and at the border with the convective zone. 
The introduction of m leads to define two classes of waves. 
Prograde (respectively retrograde) waves correspond to negative 
(respectively positive) values of m. 

Finally, all these notations allow us to define properly a criti- 
cal layer. It arises for a wave whose relative frequency crir) be- 
comes zero. The coherence with the qualitative definition given 
in the introduction is respected since cr(rc) = means that the 
excitation frequency of the wave equals -m times the angular 
velocity of the ffuid. It corresponds to a corotation resonance. 
It is important to highlight that the position of the critical layer 
depends on both wave caracteristics (m and (Tw) and shear prop- 
erties (AQ(r)). In our approximation, (rirc) = means that ^(rc) 
is constant so the critical levels are isobaric surfaces of the star. 



3. Case of the perfect fluid 

In this first approach, we treat the hydrodynamic equations as- 
suming that the fluid is neither viscous nor heat conductor. 



3.1. Equation of propagation of IGW near a critical layer 

Our aim is to calculate the Eulerian velocity field of a ffuid par- 
ticle. We introduce the time t and the unit vectors (er,ee,e(p). 
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associated with the classical spherical coordinates (r,6,(fi). The of mass conservation 
velocity field is 



y(r, 6, (pj) = r sin 6n(r)e^ + M(r, 6, cp, 0, 



(5) 



where u is the velocity associated with the wave and = 



and of energy in the adiabatic limit 



pinim = 0, 



(13) 



(14) 



rsin^Q(r) the azimuthal velocity field of the diff'erential rota- 
tion. Note that we have neglected all other large scale velocities 
such as meridionnal circulation. Then, we give the linearized 
equations of hydrodynamics governing IGWs dynamics in an The combination of these relations leads to the system presented 



- — ^ _ "I" _ br,l,m- 

P P 8 



inertial frame, i.e. the momentum, continuity and energy equa- 
tions: 

P P 

Dtp' + V.(pw) - 0, (6) 
1 p'\ 



Ti P 



—Ur = , 



where Dt = dt -\- AQd^^. The equation of energy conservation is 
obtained using the linearized equation of state 



p 



PL 
p 



r 



Tip cp' 



and assuming the ideal gaz law 

P = ^pf , 



(8) 



We define the Lagrangien displacement f as m = Dt^ and fol- 
lowing Rieutord (198 ), we expand it on the vectorial spherical 
harmonics basis (/?^, S^,T^) defined by: 



{6, if) = (6>,^)^„ 



1 



sin^ 



(9) 



' ' sm^ ' 



where are the spherical harmonics with / g M and m g [-/, /|. 
Thus: 



oo / 



(10) 



Then, we decompose p' and p' using spherical hamionics: 



oo / 

\ ^ \ ^ f /V/ 



(11) 



/=0 m=-/ 
/ 



p\r, 6,ipj) = YjTj [PiJ^^^U^^ ^)} (12) 



/=0 m=-/ 

and we obtain a new system made up of radial equations of mo- 
mentum 

PO-%-l,m = ^-^Pl^^g, 

_ Pl,m 
PO- ^H;l,m = , 

r 

^ pcr'^ir-Mm = 0, 



by Press (1981): 



dr 



dr 



= {CT^-N^)ir,l,m^ 



-yi,n 



1 d In / 2 - \ 



/(/-hi) 



P 2 

^IP 



yi,m, 

(15) 

where j/,^(r) = p[Jp. 

In stellar radiative regions, the transport of angular momentum 
is dominated by low-frequency IGWs with cr ^ N, where N is 
the Brunt- Vaisala frequency defined in Eq. (2). It allows us to 
(7) apply the anelastic approximation (Press 1981) where acoustic 
waves are filtered out. Introducing 



'i'i,m(r)=p''r^ir,i,m, 
we obtain the following equation of propagation: 



where 9i is the gaz constant and cp the specific heat per unit (j2\p 
mass at constant pressure. 



l,m 



dr2 



1(1 + 1) 



l,m 



(16) 



(17) 



1 /dlnp 
4\ dr 



1 d^ Inp 
-h - 



1 d^ln^ 



2 dr2 Ti dr2 



As the right-hand side of Eq. (17) is of order l/H^ with Hp the 
characteristic pressure or density height scale, it can be neglected 
if 



/(/-hi) 1 
^ — » 



(18) 



which is the case here. Finally, we obtain the equation of propa- 
gation of IGWs in a perfect fluid : 



l,m 



dr2 



+ 4(^)^/,m = 0. 



with 



4(^)^(^-1 



/(/ -h 1) 



(19) 



(20) 



In cartesian coordinates, this equation is called the Taylor, 
Goldstein and Synge equation (TGS). We observe that the value 
r = Tc, where cr(rc) = 0, is a singular point for this equation. 
Thus, we now focus onto the study of the behavior of this equa- 
tion around such a critical point. Then for r ^ the equation of 
propagation becomes 



^ l,m 



dr2 



/(/ + 1) Kic 



m} (r - Vc)^ 



l,m 



0, 



where 



dQ 
dr 



(21) 



(22) 
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is the value of the Richardson number at the critical level, and 

V/(/+ 1) 



(23) 



the horizontal wavenumber at the critical layer. We saw in the 
introduction (Eq. (3)) that the Richardson number is relevant to 
distinguish between the relative importance of shear and strati- 
fication. We describe this divergence of cases from a quantitave 
point of view in the following part. 

3.2. Mathematical resolution with the method of Frobenius 

The Frobenius method off'ers an infinite series solution for a 
second-order ordinary diff'erential equation of the form 

+ p(z)u + q(z)u = 0, 
where p(z) = - ^ pjz^ mid q(z) = — ^ ^jZ^ 

in the vicinity of the singular point Zo=0. For more details, a 
mathematical description of this method can be found in Teschl 
(20ii). 

In our case, p(r) = and q(r) = ^^^^^^2 - kj^c theory 
brings out two cases depending on the value of /, m and Ri^: 

- Case 1: "^-^^ Ri^ = \ 



The solution is 

r/Fro/ 



nm(^) =Mr- TcY" ^Bi(r- VcY" logdr - r,|), (24) 

where Ai and Bi are constants representing the wave 
amplitude. 

^ ./(/+!). 1 
- Case 2: ^ Ri^ 9^ - 
4 

We define the complex parameter 



1 /(/+1) ^. 

7 ^Rlc, 



and the solution is given by 



= ^2(r - r,)i/2^'^'- + B2{r - r,) 
where A2 and B2 are constants. 



\l/2-T]i^,n 



(25) 



(26) 



We exclude this special value of Ri^ and choose to consider 
only two diff'erent cases for the rest of the paper: ^^^Ri^ > \ 
and ^^^Ric < \. The results can then be extended to the case 



3.3. The case of stable critical layers: Ric > \ 

It is now time to understand the physical behaviour of the so- 
lutions given above. First of all, we propose to cut the physical 
domain into two parts: above and below the critical layer. As 
^^^Ric > |, rji^m defined by (25) is a purely imaginary number. 
In order to clearly distinguish between real and imaginary parts, 
we introduce 



0^1,, 



/(/+1) ^. 1 . 

— ^Ric - 7 = m,i 



The two solutions can be written as 



^Ffyo (r) = A+(r - r,)i/2+^-«^/.. + B+(r - VcYI^-^^^'n^ 



l,m+ 
Vi/Fro 



(r) = A_(r - rc^/^+^^i^'n + B.(r - r^yi^-i^i^^^ 



(27) 



(28) 



where ^^J^^(r) (resp. ^fj^_(r)) is available when r > Vc (resp. 
r < Vc). Booker & Bretherton (1967) and Rmgot (1*^ ) both 
explain a way to connect these solutions, considering that the 
term (r - r^)^/^^^^''" can be compared to an upward propagative 
wave of the form e^^^ and respectively that (r - r^)^/^"'^'"^ can be 
compared to a downward propagative wave of the form e~^^^ . We 
obtain the following identification : 



1/2 



upward propagating wave 



downward propagating wave 



(29) 



(30) 



In order to connect the solutions, let us observe the com- 
portement of (r - r^) above and bellow the critical level. As 
r - Vc decreases from positive to negative values, its complex 
argument changes continuously from to -n (Ringot 1998). 
Mathematicaly, we get : 



\ir>rc: {r - Vc) 



ll2±iairn — \y, _ y, \l/2±iai^m^-i7T/2^±7Tai^„ 



(31) 
(32) 



It follows that the solutions above and below the critical layer 
can be written as 



/(/+!) 



Ri^ = I without recourse to the logarithmic solution (Van y ^fj^_(^) 



/ "^IZ^^^ = - rd^/'^'""'- + B\r - r, 



Duin & Kelder 1986). 

Let us now discuss the hydrodynamical behavior corresponding 
to the situations where Ri^ > \ and Ri^ < \ . Applying 
the classical method exposed in Drazin & Reid (2004) to the gen- 
eralized spherical Taylor-Goldstein- Synge equation (Eq. 21), we 
can identify that the first regime corresponds to the case where 
the fluid stays stable with respect to the vertical shear instability 
at the critical layer. In the other case, these instability and thus 
turbulence can develop. This clearly shows how this is neces- 
sary to go beyond the current used formalisms for stellar evo- 
lution where IGWs and vertical shear instability are considered 
as uncoupled. Indeed, if a fluid becomes shear unstable, mixing 
occurs that modify the local stratification and thus IGWs propa- 
gation (see e.g. Brown & Sutherland 2007; Nault & Sutherland 
2007). Thus, we will now distinguish the case of critical layers 
when the fluid is stable from the one when it is unstable. 



l,m+ 
l,m- 



_ „ \l/2-iai^,n 



|l/2+/(2/,^ 



iBe- 



\l/2-iai^m 



(33) 

Physically, these equations can be explained this way. Starting 
from above the critical layer, the downward propagating wave 
passes through the critical layer and is attenuated by a factor 
equals to e'^^'""^. At the same time, starting from below the criti- 
cal layer, the upward propagating wave is attenuated by the same 
factor and its amplitude becomes equal to A. We also underline 
that both waves take a phase difl'erence when they cross the crit- 
ical layer. In Fig. 1, we represent the attenuation rate of difl'erent 
waves defined by the numbers / and m passing through a crit- 



ical layer. Greater is the ratio 

-4- 



1(1+1) 



stronger is the attenuation 



Att=^ " V m2 ' 4 same value of Ri^. The axis scale de- 

pends on / and m because the condition of validity of this result is 



Ri. > 



4 



. We so deduce that waves of high ratio (which 



not necessary corresponds to high order) are strongly attenuated, 
if they reach their critical layer. 
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Fig. 1. Attenuation rate Att=^ > ^ ^ of the wave passing 
through a critical layer as a function of the Richardson number. 

/(/ + 1) 

We observe that Att increases with — . 



We have not yet discussed a latter point: the choice of the method 
of resolution. Most of the publications concerning IGWs use an- 
other process to solve the equation of propagation ( ress 198 i; 
Zahn et al. 1997; Mathis 2009). In fact, the WKBJ theory is par- 
ticularly adapted to the resolution of this equation. However, it 
is not convenient in our case because it imposes a condition on 
the value of the Richardson number as demonstrated in appendix 
A. It is shown that the WKBJ approximation is available only if 
Ric » \ . Despite this restriction, let us write the solution. 
By separating the domain into two parts, we obtain 



vr/WKBJ 

Lm+ 



1 



As 



ky 



Ri » i ^ r - Tc \ 



it comes : 



m 



V/(/TT)Ri, 



-hZ)+(r - re) 



(35) 

(36) 
(37) 



It is comforting to see that both methods (Frobenius and WKBJ) 
give the same solution, to a multiplicative constant, when the 
Richardson number become high: 



rFro 
l,m 



Ri,: 



WKBJ 

l,m 



(38) 



3.4. The unstable case; Ri^ < 

In the unstable regime, the Frobenius method also gives a so- 
lution but we are not able to indentify upward and downward 
propagating waves because of shear-induced instability and tur- 
bulence. As a consequence we can not connect the solutions at 
the critical layer. In order to avoid this difficulty, we here propose 
to solve Eq. (21) following the method developped by Lindzen 



& Barker (1985). They applied it in the case of cartesian coordi- 
nates and our bringing is to generalize it to spherical coordinates. 
The parameter 7//^^ defined in Eq. (25) is real in this case and we 
introduce 

X = kHc(r-r,). (39) 



Equation (21) becomes: 

d^^/,^(X) 
dX2 ^ 



X2 



lhF/,^(X) = 0. 



(40) 



We are seeking solutions of the form ^/,^(X) 
where (^i m is the solution to the Bessel equation 



+ - inlm + ^^)^Um = 0. (41) 



_d_ 
dX 



Consequently, O/,^ is a combination of the Bessel's modified 
functions Ir]i„,iX) and I-r^^^^^X): 

^Um = ^l/,,.(^) + K2l-rj,JX). (42) 

The final solution is given by: 

^UX) = X^ {Kilrj^JX) + K2l-rj,JX)) . (43) 

We would like to calculate the reflection and transmission coef- 
ficients of the wave passing through the unstable critical layer. 
We assume that the fluid has the profile described in Fig. 2. We 
decompose the studied region into three zones defined by the 
value of the quantity ^^^Ri^. In zones I and III, the Richardson 
number is high enough to allow us to apply the WKBJ method 
described in the previous part. In unstable zone II around the 
studied critical layer however, we must use the solution with the 
modified Bessel functions. Moreover, we assume here that the 
thickness 2S of this latter is small in comparison with the char- 
acteristic length of the problem. Then, we can consider that the 
wavenumber ky is constant as 



/(/-hl)Ric /(/+!) 



6^ 



(44) 



Let us consider a wave coming from the overside of the critical 
layer (zone I). It is partly transmitted toward zone III and partly 
reflected backward zone I. We so write the solutions correspond- 
ing to the three zones: 



^-iky{r-rc) _|_ j^^ikv(r-rc) ^ 



^iii(r) = re*"*'-'''', 

(45) 

where is available if (r - rc) > S, ^ii(r) if -6 < (r - rc) < 6 
and ^in(r) if {r - rc) > -6. The coefficients A, B, T and R are 
calculated thanks to the following four continuity equations: 



df/ 



T„(r, + 5), 

{rc + 6)=^{rc + 6), 
ar ar 
'i'luirc -S) = >P//(r, - S), 



(46) 



dr 



dr 



They correspond to the continuity of the solution and of its 
first derivative, which physically means that both displacement 
and mechanical stresses are continuous. After some algebra. 
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^ -r=rc+5 



Fig. 2. Assumed neighbourhood of an unstable critical layer for 
the calculation of the IGWs' reflection and transmission coef- 
ficients. We assume that the unstable region around the criti- 
cal layer (in red) has a thickness given by 26 (Zone II). The 
surrounding regions where IGWs are propagative are in beige 
(Zones I and III). 

we obtain the coefficients R and T of reflection and transmis- 
sion, depending on the stiffness 6, on the vertical and horizontal 
wavenumbers kvc and Uhc (see Eq. (44) and (23)) , and on the 
variable (Eq. (25)). In order to lighten the formula, we note 



instead of I+t^i^X^hc^)- Then, we obtain 



R = 



Rn 



^denoml "I" ^denom2 



(47) 



(zone II), arbitrarily choosen, that points the non-local charac- 
ter of unstable turbulent layers. Let us underline the main re- 
sult: both coefficients are greater than 1 when ^^^Ri^ is small 
enough. Consequently, for a low Richardson number at the criti- 
cal layer, the wave can be over-reflected and over-transmitted at 
the same time. It means that on the contrary of the first stable 
case, the wave take potential energy from the unstable fluid and 
convert it into kinetic energy. In other words, the turbulent layer 
acts as an excitation region. If \R\ < 1 and \T\ < 1, we speak 
about IGWs "tunneling" (Sutherland & Yewchuk 2004; Brown 
& Sutherland 2007; Nault & Sutherland 2007). Another remark 
concerns de dependency of \R\ and |r| with m. We denote that 
the sign of m does not matter since only its square appears in 
the expressions. Physically, it shows that the critical layer's ac- 
tion is the same on prograde and retrograde waves. This point is 
of importance because we know that other dissipative processes 
occuring during the propagation of IGW discriminate between 
both types of waves. In order to visualize the action of the criti- 
cal layer on differents waves. Fig. 2 shows the level lines |^| = 1 
for 1 < / < 5 and 1 < m < /. We previously said that it is useless 
to consider negative values of m since \R\ depends only on m^. 
The pair (/,m) is indicated on each ligne, followed by the value 
of Lines in the same color correspond to the same value 
of /. These lines mark out the limit to observe over-reflection, 
for a chosen wave. We observe that higher is the value of 
stronger is the condition on Ri^ to observe an over-reflection. 

3.5. Choice of the method 

We have decided to apply different methods to solve the stable 
and unstable cases. However, it could be legitimate to wonder if 
both methods are equivalent from a mathematical point of view. 
In this part, we present a short comparison between the solu- 
tions obtained with the method of Frobenius and the one with 
Bessel functions. The modified Bessel function Irji^(X) can be 
computed using 



with 



kHcl'rj,,n-^\^-ikvc]lr^i,„, 



R. 



denoml - '^Hc^-rii,m^-^lm 
2kvc 



-^denom2 

and 
with 



cos (7//,^7r) -h ^ {lmJ-m,m + ^-mJmJ 



^denom 



(48) 



_ 2ikvc 

^ num ~ c ■> 

on 

Td enom — ^denoml "I" ^denom2- 

We have now calculated the transmission and reffection coeflS- 
cients of an IGW, under some assumption, through an unstable 



region around a given critical layer where Ri^ < 



1 

4/(/+l) 



. We rep- 



resent the level lines of \R\ and |r| in Fig. 3 for / = 4 and m = ±3. 
They are plotted as function of the Richardson number at the crit- 
ical layer Ri^, growing from to its maximum value defined by 

0.11. 



^Ricm«x = i - For instance, in Fig. 3, Ri, 



~ 4 4(4+1) 

The other variable is the half- thickness S of the unstable layer 



(1 \m,m ^ 
^ ' k=0 



(49) 



At the neighbourhood of the critical layer, X tends to and the 
first-order expression is 



X 



1 



r(^/,m + 1) 



(50) 



Close to the critical layer, the global solution given in Eq. (43) is 
then 



rBessel^^^ = | Ai — — -h ^ + 0(X^) . (51) 



Lm 



Now, let us remind the expression of the solution given by the 
method of Frobenius (Eq. (26)), rewritten with the previous no- 
tations 

In conclusion, for a fixed couple (/,m), 'i'f^ir) and ^^f^'^XX) 
vary in the same way as function of r - r^. 
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Fig. 3. Level lines of reflection \R\ and transmission |r| coefficients of an IGW at a critical layer as a function of the Richardson 
number Ri^ and of the thickness S of the critical layer (zone II in Fig. 2). The top panels represent the level lines for an arbitrary 
choosen value of (/, m) = (4, ±3), while the bottom panels show levels lines \R\ = 1 and |r| = 1 for diff'erent couples (/, m). 



4. Case of the non-perfect fluid 

From now, we have studied the role of the critical layers as- 
suming that the fluid was perfect. In order to make the problem 
more realistic we include in the second part the viscosity v of the 
fluid and the coeflficient of thermal conductivity k (e.g. Koppel 
1964; Hazel 1967; Baldwin & Roberts 1970; Van Duin & Kelder 
1986). 

4.1. Equation of propagation of IGW near a critical layer 
The linearized equations of hydrodynamics in Eq. (6) become: 



DtU 



Dtp' + V.Oom) = 0, 



(52) 



The following method for the building of the propagation equa- 
tion is adapted from the work of Baldwin & Roberts (1970). We 
assume that the mean density p of the fluid nearly takes a con- 
stant value, that is to say that | ^ is small compared with others 
characteristic lenghts. First, we project Eq. (52) onto Cr and ap- 
ply the operator V^. Then we apply Dt - kV^. After combination 



with the two other equations, we obtain : 



(A-^V^) 



(A - v^^)(y^Ur + (Q!' + -0!)d^Ur 



(53) 

As done in the previous part, we decompose the radial velocity 
on the basis of spherical harmonics 



l,m 



(54) 



where cr = cr^ -\- mAQ(r). Moreover, it is easier to work with 
dimensionless numbers. For this reason, we introduce the nota- 
tions detailed in Tab. 1 . 



Table 1. Dimensionless numbers used for the resolution of Eq. 
(53). L and V are respectively the length and velocity scales. 



Pr 


v/k 


Prandtl number 




VL/y 


Reynolds number 


Ri 


(LN/Vf 


Richardson number 
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Then, Eq. 53 becomes for each pair (/,m) such as / g M and 
m G [-/, /I : 

(A/ -kl- icrR^Pr) (A/ -kjj- io-R^) (A/ - k]j) Ur,l,m 

= -klRlP,NX,i,^, (55) 
where A/ is the scalar spherical Laplacian operator : 

Ai = dl+-dr-^-^^. (56) 

Lastly, we introduce r] = {imOf^R^Y'^ir - r^) and thanks to a 
developpement close to the critical layer we obtain the sixth- 
order equation: 




1 d%,^ 



+ n 



/(/+!) 



1 \ d^Xi,m 2 d^xi,r. 



dif- 



PrI ^7/4 P, drf 



Fig. 4. Left: curves (C\,C2,C^) defining the basis 
{UuU2,U^,VuV2.y?>). Right: curves {C[,C^,C^) defining 
the basis (wi, ^2, W3, vi, V2, V3). 



(57) 



Thanks to this solution, we obtain the expression of the radial 
where no^ Xi,m = P^'^nKum- Equation (57) can be compared displacement of the wave: 
with the one obtained by Press (1981) who has neglected the 
viscosity v. Moreover, we can see that if we take v = /c = 0, and 
without forgetting that Py depends on v too, Eq. (57) is identical 
to Eq. (21) for the perfect fluid. 



(63) 



4.2. Mathematical resolution 



we will clarify the function /(k^Py) of the thermal difl'usivity 
The resolution of Eq. (57) requests several substitutions and coefl&cient later. 

is quite complex. The detailed calculation can be found in ap- There is still the last step to get over. We must find a curve along 
pendix B and we give here only the main steps. We draw our in- which we integrate the solution, that is to say, we determine a 
spiration from Hazel (1967); Baldwin & Roberts (1970); Koppel and b. During the calculation detailed in the appendix B, Eq. 
(1964); Van Duin & Kelder (1986) who solve the same equation (104): 
in cartesian coordinates. However, the resolution in spherical co- 
ordinates has some diff'erences. The aim is here to rewrite the 
equation under a known form: the Whittaker diff'erential equa- 
tion; then, after some algebra, we can show that Eq. (57) can be 
written in the following form 



1 -h 



1 



dt 







(64) 



has not been used yet. A sufficient condition to make this egality 



\-^lm A 
2 + - 



^ 00. Considering that r 



y = 0. 



(58) 



true IS 

complex argument of t), we have 



fe^'^' (6t being the 



3 \t\- 



00 <^ Of = 



2n 



(65) 



The solutions of Eq. (58) are thus the Whittaker functions 
(Abramowitz & Stegun 1965): 

There are several possibilities for the choice of the curves. 

1 _ _ \ Koppel (1964) propose to build a basis of six solutions using 

the curves (Ci,C2,C3) represented in Fig. 6 (left). Therefore, the 
solutions of the sixth order equation (Eq. (57)) are linear combi- 
where is the confluent hypergeometric function of Kummer nations of 



1^1 



+ M/,^+A;l+2M/,^;^ , (59) 



iFi(a;b;z) = ^ 



(b)nn\ 



(60) 



Ut(r]) = f(K,P,) 

JCi 



^V3^3M/,^-3/2 



lFi\-Mi^rn'A-^^^l,mU]dt 



(66) 



n-l 

with(k)n = n(^+o 

Gamma function. 



T(k+n) 

nk) 



and (k)o = 1, r (z) being the usual 



and Vi(ri) with / e {1,2,3}, where Viirj) corresponds to Ui(ri) 
with the opposite sign for Mi^m- 



2+^^Ri, 



(61) 



and 



4.3. Application to IGWs 



It is now time to apply these mathematical results to IGWs. With 
this aim in view, let us remind the solution for the perfect fluid 
obtained in the first part. For the moment, it is not necessary 
(62) to distinguish between the stable regime, i.e. ^^^Ri^ > t, and 



8 



L. Alvan, S. Mathis, T. Decressin: Critical layers for internal waves in stellar radiation zones 



the unstable one, i.e. ^^^Rir < \. The Frobenius solution is 



1(1+1) 



Rir as a corn- 



available everywhere if we take rji^m = -yjl ~ 

plex number. Remembering that = P^r^h\i,m{r), the ra- 

dial Lagrangian displacement is 



(67) 

In order to make the reading easier, we have removed the in- 
dices r; /, m. Indice P designates the solution for the perfect fluid. 

Concerning the solution for non-perfect fluids. Hazel (19b / ) pro- 
poses to use another basis of the solutions of Eq. (57). The 
curves {P[,C2,C'^) are represented in Fig. 6. The new basis 
(^1,^2, ^3, vi, V2, V3) is defined by: 

- ui = Ui + U2 + U^ and vi = Vi -h ¥2 + Vs, 



- U2 

- U3 



U1 + U2 + u; and V2 = Vi + y2 + v;, 
Vi + y* + ^3, 



Ui-\- U*-\- U3 and V3 



where X* is the complex conjugate of X. Consequently, the so- 
lution for the non-perfect fluid (indice NP) can be written as : 

1 ^ 

^NP± = —TIT 

^ ^ i-l 

A relation between both solutions |p+ and |np± exists if we con- 
sider that the Reynolds number is great. In the case of the Sun 
(e.g. Brun & Zahn 2006), the microscopic viscosity in the radia- 
tive zone is weak. For this reason, it is appropriate to consider 
that the Reynolds number Rq = ^ is much greater than 1 . This 
assuption leads to the relation 



^NP± „ — ^ ^p+- 



(69) 



Baldwin & Roberts (1970) give tables for the asymptotic com- 
portement of Uf and v/ ( ie 1,2,3). The solutions ui, U2, vi et 
V2 diverge when Rq -\-oo. They are therefore physically unac- 
ceptable and we deduce that 



1 



,NP± - 



(70) 



Tables (2) and (3) give the asymptotic expressions of U3 and V3 
as a function of the sign of r - r^. 



In the stable case where ^^^Ri^ 



> 4. 



J^Ri, 



is a real number. The same findings than in the previous part 
can be made: after the passage through the critical layer, the 
wave is attenuated by a factor e'^^'"^^. In the unstable case where 

^^^Ric < |, ai^m = -yj^-^^Ric - \ is a complex number, so we 
can not interpret the solution as upward and downward propagat- 
ing waves. But the expressions given in Tab. (2) and (3) remain 
comparable to those in the first part and we deduce that the cal- 
culation of R and T will lead to the same result: the possibility 
of over-reflection and over-transmission. 

4.4. Radiative and viscous dampings 
4.4.1. General equations 

We volontary left aside the factor /(/c, Pj-) in Eq. (66). In order 
to establish its expression, Zahn et al. (1997) use the equation of 



Table 2. Expressions of U3 and V3 when Rq -hoo, below the 
critical layer. 



W3(^) 



Table 3. Expressions of U3 and V3 when Rq -hoo, above the 
critical layer. 



r <r. 



V3(^) 



2n 



In 



r(f + iairn) 



IGWs propagation taking into account heat difl'usion but with a 
viscosity coefficient (v) equals to zero. Here, we generalize their 
result for y 9^ (i.e. ^ 0) and we obtain 

-T(/C,Pr)/2 



f(K,P,) = e- 



where 



T(^,Pr) =[/(/+!)] 



r 



Kd+Pr) 



(71) 



(72) 



We have introduced the general expression for N^, the Briint- 
Vaisala frequency, to be able to take into account chemical strat- 
ification. Then, we have 



with A^^ 



(73) 

jf, (^ad - V) and N^ = jf^V^ where Hp = |dr/dlnP| 

is the pressure height-scale, V = (dlnT /dlnP^ the temperature 

gradient and = (51ny[Z/51nP) the mean molecular weight (jj) 
gradient. Moreover, we have introduced the generalized equation 
of state (EOS) given in Kippenhahn & Weigert (1990): 



dp _ I dP _~&T ^ ^dfi 



(74) 



where 6 = - (d\np/dlnT)p^ and = (^lnp/(91nyu)p 7^. Next, 
(T = (Tw + mAQ is the Doppler- shifted frequency of the wave 
relative to the fluid rotation with an excitation frequency cr^. Vc 
and rcz are respectively the positions of the critical layer and 
of the boundary between the studied radiative zone and the con- 
vection region where IGWs are initially excited. This damping is 
another source of attenuation independent from the presence of a 
critical layer. Moreover, as shown in §3.4. and §4., we will have 
to consider IGWs reflected and transmitted by unstable critical 
layers in addition to those initially excited by convection. Then, 
we introduce the following notation 

T[/c,Pr,ri,r2] = 



[/(/+!)] 



rf 



cr4 lAf2-cr2/ r 



1 



dr. 



(75) 
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where ri and r2 are respectively the studied IGW's emission 
point and the studied position with ri > r2 (in the opposite case 
where ri < r2, limits in the integral have to be reversed). This 
will enable us to describe radiative and viscous dampings in ev- 
ery configuration. Note that in stellar radiation zones ^ I 
(Brun & Zahn 2006) inducing that the damping is mostly radia- 
tive. We have now to compare it with critical layers' eff'ects. 

4.4.2. Progrades and retrogrades waves 

For a same environnement, prograde waves have a frequency 
lower than retrograde ones (e.g. Eq. 4). We choose for exam- 
ple a couple of IGWs with the same excitation frequency, ctq, 
same number / and opposite azimuthal degree m. Thus, we com- 
pare a prograde wave of frequency (Ti(r) = ctq - |m|AQ and 
a retrograde one of frequency (T2(r) = ctq + |m|AQ. We obtain 
(Ti(r) - (T2(r) = -2\m\ (^(r) - Qcz) < in the presence of neg- 
ative gradient of and (Ti(r) - cr2(r) > if the gradient is posi- 
tive. Now, let us underline that r given in Eq. (71) varies globaly 
as ^ . As a consequence, assuming that a negative Q-gradient is 
present near the excitation layer (see §. 6), the radiative damping 
is stronger for prograde waves than retrograde waves. Therefore, 
prograde IGWs are absorbed by the fluid much closer to their 
region of excitation while the retrograde waves are damped in a 
deeper region. This process is responsible for the net transport 
of angular momentum by IGWs in stars. On the other hand, crit- 
ical layers doesnt introduce such net bias between prograde and 
retrgrade IGWs because their efl'ects scale with m^. 



4.4.3. Dependency in / and m 

The second remark concerns the variation of radiative and vis- 
cous dampings as a function of / and m. On one hand, looking 
at the multiplicative factor, which is in front of the integral in 

Eq. (72), we can roughly write that r oc ^^^^^ — . On the other 
hand, the expression of the attenuation factor due to stable crit- 

(1(1+1) \^/^ 
-^^^1 , which is always greater 

than unity since \m\ < I. The comparison between radiative and 
viscous dampings and the one due to stable critical layers for an 
IGW with given (/, m) will be examined in §6.2.3. 

4.4.4. Location 

Lastly, radiative and viscous dampings occur throughout the 
whole propagation of IGWs. On the contrary, critical layers are 
localised. Moreover, there is a condition for a wave to reach a 
critical layer: the rotation speed of the fluid must be of the same 
order than the wave frequency to hope observing cr = 0. Finally, 
all IGWs are concerned by radiative and viscous dampings, 
which increase around critical layers since r oc cr"^. 

In Tab. (4), we sum up the diff'erent cases studied in these work. 
Depending on both properties of the fluid and of the studied 
wave, these latter is submitted to diflferent phenomena. 



5. Transport of angular momentum 

As emphasized in the introduction, our goal is to study the trans- 
port of angular momentum in stellar radiation zones and to un- 
ravel the role of critical layers. Therefore, the first step is to recall 
the flux of angular momentum transported by propagative IGWs 
and by the shear-induced turbulence. Moreover, to illustrate our 



stable 
Critical Layer 




Initial excitation region 
Shear unstable critical layer 



Damping (Th.) and 
possible over-reflection 
& over-transmission (CL) 



Initial wave kinetic energy flux 

Reflected or/and transmitted wave kinetic energy flux 



Fig. 5. The two studied configurations in a low-mass star with an 
external convective envelope. Left: the case of a stable critical 
layer (CL) where IGWs are damped. Right: the case of an unsta- 
ble critical layer where IGWs can be over-reflected/transmitted. 



purpose, we will here focus on the case of a low-mass star where 
IGWs are initially excited at the border of the upper convective 
envelope (see Fig. 5) by turbulent convection (Garcia Lopez & 
Spruit 1991; Dintrans et al. 2005; Rogers & Glatzmaier 2005; 
Belkacem et al. 2009; Brun et al. 2011; Lecoanet & Quataert 
2012) and by tides if there is a close companion (Zahn 1975; 
Ogilvie & Lin 2007) with an amplitude A (corresponding results 
for massive stars with an internal convective core can be easily 
deduced by reversing signs). 



5.1. Angular momentum fluxes 

5.1 .1 . Angular momentum flux transported by propagative 
IGWs 

First, we have to calculate the angular momentum flux trans- 
ported by a propagative monochromatic wave over a spherical 
surface. It is given by the horizontal average of the Reynolds 
stresses associated to the wave (e.g. Zahn et al. 1997): 



'Fj-l,m,cr (r) = pr sin ^ < Ur,l,mU^-l,m >e,^ , 



(76) 



where < ... 
shows that 



>. = -L r ...smedOd(f. Besides, Mathis (2009) 
47r 



-m 



J;l,m,cr 



(77) 



where TE;i,m,cr is the horizontal average of the energy flux in the 
vertical direction expressed by Lighthill (1986) as 



So finally, the angular momentum flux is given by 



(78) 



(79) 



To calculate this angular momentum flux, we need the expres- 
sions of Ur-i,m and p'^^. Ur-i,m IS immediately accessible from Eq. 
(33), solution of Eq. (21). It is a little more complicated for p'^^ 
because we need to go back in the calculation leading to Eq. 
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Table 4. Summary of the different studied cases. The observed effects depend on the properties of the wave (the frequency (cr) and 
degrees (/, m)) and of the fluid (thermal conductivity (/c), viscosity coefficient (v), Richardson number (Ri)). 



wave and fluid properties 




o-(rc) 


K 


y 


^Ri, 






^0 


not a critical layer 













> 1/4 


attenuation 


Fig. 1 


< 1/4 


possible over-reflection + possible over- transmission 


Fig. 3 


^0 


« 1 


> 1/4 


attenuation + radiative damping 




^ 1 


< 1/4 


possible over-reflection + possible over-transmission + radiative damping 


» 1 


non stellar case 



(21). The expression of p'^^ results from the first equation of the 
system given in Eq. (15), which can be reduced to 



d 
d^ 



|^j = (cr^-A^2)|,,,^, (80) 



applying the anelastic approximation and neglecting terms of 
order 1/L^. 

The total vertical flux of angular momentum transported by 
propagative IGWs is then given by 



rj(r)= Yj'^J-Mm,o 



(81) 



l,m,cr 



and we define the associated so-called action of angular momen- 
tum 

£j(r) = Anr^Tj. (82) 



5.1 .2. Angular momentum flux transported by shear-induced 
turbulence 

In the case of shear-unstable regions, such as region II in the 
regime where Ri^ > IGWs are unstable (e.g. jJrazm & 

Reid 2004). After the non-linear saturation of the instability, a 
steady turbulent state is reached. Then, as established for exam- 
ple in Zahn (1992); Talon & Zahn (1997), and now confirmed by 
numerical simulations (see Prat & Lignieres 2013), the vertical 
flux of angular momentum transported by shear-induced turbu- 
lence is given by 



5.1 .3. Equation of transport of angular momentum 

Let us now refocus these results in the wider frame of the com- 
plete angular momentum transport theory. Considering the other 
transport mechanisms we presented in introduction, the angu- 
lar momentum transport equation taking into account merid- 
ional flows, shear-induced turbulence and IGWs (e.g. Talon & 
Charbonnel 2005; Mathis 2009) becomes 



dr-Lj (r) where IGWs are propagative 



-\-< 

for unstable regions. 



(85) 



rT;v(r)=pr^yv(r)dr^, 



(83) 



The first term in the right hand side corresponds to the advection 
of angular momentum by the meridional circulation, where 
Ur = U2P2(cosO) is its radial component. ^ = ^ + is 
the Lagrangian derivative that takes into account the radial 
contractions and dilatations of the star during its evolution, 
which are described by rir. According to our hypothesis, we 
do not take into account the transport by the Lorentz force, 
associated with magnetic fields in stellar radiative zones (Mathis 
& Zahn 2005). 

The major diff'erence with previously published equations is that 
IGWs and shear-induced turbulence transports of angular mo- 
mentum are not summed linearly since they are, as we demon- 
strated before, intrisically coupled. Therefore, for stable regions, 
one must take into account IGWs' Reynolds stresses (Eq. 82) 
only, while for unstable regions, only the vertical turbulent flux 
given in Eq. (90) must be taken into account. Let us now exam- 
ine each case of critical layer. 



where 



Ri^ 



NII{k + Vh) + NIIvh 



(84) 



Ri^ = 1/6 is the adopted value for the critical Richardson num- 
ber and Vh is the horizontal turbulent viscosity for which we 
assume the prescription derived by Zahn (1992). 



5.2. Stable critical layer \Ric > 
5.2.1 . Case of the perfect fluid 



1 

4/(/+ 1) 



This is the simplest case of a stable critical layer in a perfect 
fluid. Then, we apply Eq. (79) to obtain the expressions for the 
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transported fluxes by propagative IGWs below and above the 
critical layer: 



(86) 



where A is the initial amplitude of the IGW at r = rcz, cti^, 



^'■^Ric - I (see Eq. (27)) and 



Ji,, 



--< fPf(C0S(?)f >fl= 

2 (/+M)! 



sin^d6> 



2/+l(/-N)!' 



(87) 



where PJ^ are the associated Legendre polynomials. Booker & 
Bretherton (1967) have obtained similar results in cartesian coor- 
dinates. Let us make two remarks about these expressions. First, 
we see the expected attenuation of the flux by a factor e~^^^'^"' 
when the wave passes through the critical layer. Moreover, 
'J^j\i,m,cr± depends on m (and not on m^). As a consequence, we 
recover the classical result that prograde waves (m < 0) and 
retrograde ones (m > 0) have opposite angular momentum flux 
(respectively a deposit and an extraction). Finally, the monochro- 
matic action of angular momentum X/;/,m,o-(^) = 47rr^^F/;/,m,o- is 
constant in each region because of the absence of dissipation. 

5.2.2. Case of the non-perfect fuid 

We saw in the previous part that the solution of the equation of 
propagation in the case of a non-perfect fluid is similar to the 
one obtained for a perfect fluid. For this reason, we are allowed 
to apply the same method for the calculation of Tj-im,cr and we 
obtain 



-T [k, Pr, rcz, r] 



Tj-i,m,cT (r < rc) = 
1 1 



(88) 



J"/,m ^-27Tai ^^-T [k, Pr, Tcz, r] 



The difference with Eq. (86) comes from the introduction of 
radiative and viscous dampings (Eq. 75). 



(§5.2.2.), we obtain: 

Tj-i,m,cr {rc + d<r< rcz) = 

L ^rnA^ _3um_ -T [k, P,, rcz, r] 
r^2 /(/+1) 
1 1 

r^ 2 



xe 



^2|^|2^-T [k, Pr, rcz, rc + S] 
T [k, Py, r, rc + S] 



>Jl,n 



/(/-hi) 



(89) 



where we identify the transport induced by the incident wave, 
which propagates downward, and the one induce by the reflected 
one, which propagates upward. Then, we can see that the unsta- 
ble critical layer is a second excitation source for IGWs propa- 
gating in region I and that the angular momentum transport will 
be modified in such situation, particularly when over-reflection 
(\R\ > 1) occurs. 

5.3.2. Region \\ \ rc - 6 < r < rc + 6 

In this unstable region, we directly use results obtained in §5.1.2. 
to describe the flux of angular momentum transported by the 
shear-induced turbulence, i.e.: 



T^T-y (rc - S < r < rc -\- S) = pr^vydrfl. 



(90) 



5.3.3. Region III: r>rc-S 



Using the same methodology that for region I, we get in a 
straightfoward way: 



Tj-i,m,cT (r>rc-S) 



1 1 

2 



xe 



T [/c, Pr, rc - S, r] 



/(/-hi) 



(91) 



where we indentify the transport induced by the transmitted 
wave, which propagate downward. As in region I, we can 
see that the unstable critical layer constitutes a secondary 
excitation source for IGWs propagating in region III and that 
the angular momentum transport will be modified, particularly 
when over-transmission (|r| > 1) occurs. 

Since the general theoretical framework has been given, we have 
now to explore the possibility of the existence of the two differ- 
ent regimes (stable and unstable) along the evolution of a given 
star. As a first application, we choose to study the case of a solar- 
type star which has already been studied without critical layers 
by lalon & Charbonnel (2005). 



The conclusion is that, in the stable case, the attenuation due to 
the passage through a critical layer is added to dampings due 
to dissipation. This observation will lead us to simply imple- 
ment the role of stable critical layers as an additional term in the 
damping coefficient. 



5.3. Unstable critical layer \Ric < - 




5.3.1 . Region I: + ^ < r < rcz 

Using the summary given in Fig. 5 for the unstable case and 
the obtained results for regions where IGWs are propagative 



6. A first application: the evolution of a solar-type 
star 

6.1. The STAREVOL code 

We use the one dimensional hydrodynamical Lagrangian stellar 
evolution code STAREVOL (V3.10), and the reader is referred 
to Lagarde et al. (2012) and references therein for a detailed 
description of the input physics. We simply recall the main 
characteristics and parameters used for the modelling that are 
directly relevant for the present work. We use the Scharwschild 
criterion to determine convective zones position, and their 
temperature gradient is computed according to the MLT with 
a o^MLT = 1.75. The solar composition is taken from Asplund 
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et al. (2005) with the Ne abundance from (Cunha et al. 2006). 
We have generated the opacity table for temperature higher 
than 8000 K following Iglesias & Rogers (1996) by using their 
website^. Opacity table at lower temperature follows Ferguson 
et al. (2005)^. The mass loss rate is determined following 
Reimers (1975) with a parameter tjr = 0.5. The increase of mass 
loss due to rotation is taken into account following Maeder & 
Meynet (200 ). However due to the small mass loss and velocity 
of our model this effect remains weak. 

In radiative regions, we follow Mathis & Zahn (2004) formalism 
for the transport of angular momentum and of chemicals as well 
as the prescription from Talon & Zahn (1997) for the vertical 
turbulent transport (Eq. 84). We assume that convective regions 
rotate as solid-body. The treatment of IGW follows laion & 
Charbonnel (2005, 2008) with the difference that the volumetric 
excitation by Reynolds stresses in the bulk of convective zones 
(e.g. Goldreich & Kumar 1990; Goldreich et al 1994; Belkacem 
et al. 2009) is consistently computed at each time-step as a 
function of their physical properties. 

We start with an initial model of 1.0 M© at solar metallicity, 
with an initial surface velocity of 70 km s"^ The rotation profile 
is initially flat. We add magnetic braking through the following 
law: ^ = -KQ"^ with a constant K = 3x 10^^. This value has 
been calibrated to reproduce the surface velocity determined in 
the Hyades by (Gaige 1993). 

6.2. The effects of critical layers 
6.2.1 . Location of critical layers 

We theoretically studied the impact of the critical layer for 
a given IGW, thus assuming that there are some radii where 
the relation cr(r) = is satisfied. As a consequence, the first 
question we may answer thanks to the simulation concerns 
the existence of such critical layers and their location in the 
radiative zone. Figure 6 shows that critical layers do exist in 
the studied solar-like star's radiative core. In the three panels, 
we superimpose the rotational velocity of the star's interior 
as a function of the normalized radius with the position of 
potential critical layers, marked out with colorized squares 
which correspond to positions where cr = (Tw + mAQ = 0. Each 
panel corresponds to a given value of the excitation frequency 
cr^. As expected, the positions of the critical layers only depend 
on the azimuthal number m of the wave, and not on the degree 
/. Thanks to this plot, we confirm that critical layers exist in the 
studied solar-like star. However, we have not already taken into 
account the fact that all waves cannot reach these positions. 

First of all, we only take into account retrograde waves because 
the prograde ones are damped immediately after their excitation 
(see § 4.4.2). Moreover, in the left panel of Fig. 7, we have rep- 
resented in logarithmic scale the luminosity of a given IGW at 
the moment of its initial stochastic excitation by the turbulent 
convection as a function of its degrees / and m, following the 
spectrum adopted in Talon & Charbonnel (2005). The difference 
between the three plots is the value of the excitation frequency 
(Tw. Thanks to this representation, we can see that the maximum 
of excitation lies in a domain where / and m are close and quite 
small. Moreover, it shows that for each excitation frequency, the 



^ http : //adg . llnl . gov/Research/OPAL/opal . html 
^ http : //webs . wichita . edu/physics/opacity/ 
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Fig. 6. Position of the critical layers for retrograde IGWs for 
three different frequencies, superimposed with the rotational ve- 
locity profile in the star. 

amplitude of the excited wave depends on / and m, and may be 
close to zero. As a consequence, some critical layers represented 
in Fig. 6 may belong to a non-excited wave, or to a wave com- 
pletely damped at this depth. Fortunately, obtained results show 
that some waves really meet their critical layer. 

6.2.2. Interaction between waves and critical layers 

Concerning the way these latters interact with the surrounding 
fluid, the theory predicts two possible regimes depending on the 
value of the Richardson number Ri^ at the critical level. In our 
simulation, it appears that for every detected critical layer, the 
relation Ri^ > \ , which correponds to the stable regime, is 
verified. As a consequence, we establish that the second unstable 
regime (with the associated possible tunneling or over-reflection 
and transmission) does not occur for the solar-like star of our 
simulation; forthcoming studies may explore other types of stars 
at different evolutionary stages, to see if this regime can occur. 
Therefore, for the solar- type star studied here, we only imple- 
ment in STARE VOL the terms related to stable critical layers: 
each time a wave passes through a critical layer, it is thus damped 

with a coefficient e > ' ^ (see Eq. 88). 

6.2.3. Effect of critical layers 

Since all interactions between waves and critical layers are of the 
same kind in the studied star, we can concentrate on the quanti- 
tative importance of their effect on the transport of angular mo- 
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Fig. 7. Left: Luminosity of the waves at the location of their excitation as function of / and m. Right: Ratio between tcl and 
Ttot = Tcl + Trad (scc Eq. (92)). 



mentum. We know that in the stable regime (§5.2.2.), the wave 
passing through its critical layer is damped by the factor given 
in Eq. 88) that is added to radiative damping (here ^ I and 
the viscous damping is thus negligible), which has already been 
taken into account in previous works (e.g. Talon & Charbonnel 
2005). In the right panel of Figure 7, we thus choose to represent 
the ratio between tcl, the rate of attenuation due to the passage 
of the wave through its stable critical layer and the sum Ttot of 
this rate and the one of the radiative damping. The explicite for- 
mula is 



Tcl 

Ttot 



(92) 



27T 



/ /(/+!) 
V m2 



Ri^ 



27r 



V m2 



Ri^ 



+ [/(/+ 1)]^ 



nrzc 

J ' 



1 



rdr 



On the contrary of the left panel, here are represented only waves 
which meet a critical layer. That is the reason why several white 
zones are seen. They correspond to the waves which have been 
attenuated before reaching the depth of their critical layer. In the 
case of low (Xw (upper panels in Fig. 7), the high degree waves 
(/ > 45) are simply not excited, as we can see on the right. In red 
zones, the role of critical layers is important in comparison with 
the radiative damping while dark blue regions are those where 
the radiative damping dominates . Therefore, this figure shows 
that critical layers should be taken into account. 



6.2.4. Evolution of the rotation profile 

Let us now concentrate on the evolution of the rotation profile 
when momentum deposition due to the interaction between 
waves and critical layers is taken into account. The other 
transport mechanism (IGWs' radiative damping eff'ects, the 
shear-induced turbulence and the meridional circulation) have 
been previously implemented in the code ( Talon & Charbonnel 



In the center of the star, the rotation velocity is lower and 
increases the radiative damping of retrograde waves (see Eq. 
72). This forms an angular momentum extraction front which 
propagates from the core to the surface to damp the diff'erential 
rotation. Three fronts are seen in the top panel of Fig. 8 and 
have the same form that those already obtained in Talon & 
Charbonnel (2005). To isolate the action of critical layers on the 
evolution of the rotation profile, we have superimposed in the 
bottom panel of Fig. 8 the surface velocity as function of the 
evolution time with (black line) and without (red line) critical 
layer eff'ects between 2.8 x 10^ and 3.5 x 10^ years. Both curves 
are nearly identical wich shows that despite their local action 
showed in Fig. 6, IGWs' stable critical layers do not disturb 
the dynamical evolution of surface velocity in the case of the 
studied star. This can be easily understood since the radiative 
damping becomes mostly efficient around critical layers' posi- 
tions because of its dependance on cr~^ . Moreover, if in Talon 
± Charbonnel (2005), the IGWs action and the shear-induced 
turbulence have been added as uncoupled physical mechanisms, 
it has been demonstrated that in the upper region (r > r^) where 
IGWs are propagative, the coefficient vy is negligible, while in 
the inner one (r < r^) the diff'erential rotation has been damped, 
leading to the same result with a transport dominated by the 
meridional circulation. 



This clearly indicates that unstable critical layers will lead to ma- 
jor modification of the transport of angular momentum in stel- 
lar interiors. A systematic exploration of diff'erent type of stars 
for diff'erent evolutionary stages will be undertaken in a near fu- 
ture to explore their possible existence. Moreover, in order to 
give quantitative informations, it is necessary to improve the way 
waves are excited in this model. A future work will implement 
a prescription about penetrative convection processes. The best 
way to do it is to use a realistic numerical simulation of such 
mechanism to obtain the excitation spectrum at the base of the 
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Fig. 8. Evolution of the rotation profile where the role of the crit- 
ical layers in the transport of angular momentum is taken into 
account. The curves are labelled according to the corresponding 
ages in Gyr. Parameters are indicated above where K is the brak- 
ing constant. Right : Comparison between the evolution of the 
rotation with (red) and without (black) taking into account the 
eff'ect of the critical layers. 

convective zone (Rogers & Glatzmaier 2005; Brun et al. 2011; 
Alvan et al. 2012). 

7. Conclusion 

In this paper, we study in details a new mechanism of inter- 
action between IGWs and the shear of the mean flow that oc- 
curs at co-rotation layers in stably stratified stellar radiation 
zones. Taking advantage of the work realized in the literature 
about atmospheric and oceanic fluids, we highlight the similari- 
ties with such stellar regions and propose an analytical approach 
adapted to the related case of deep spherical shells. Then, the 



use of spherical coordinates brings difl'erences in the equations 
and make their resolution more complicated but the final results 
are comparable. We then demonstrate the intrinsic couplings be- 
tween IGWs and the shear-induced instabilities and turbulence 
that can thus not be added linearly as done previously in stel- 
lar evolution literature. Then, we highlight the existence of two 
regimes where the interactions between IGWs and the shear at 
critical layers are strongly diff'erent: 



- in the first case, the fluid is stable and IGWs amplitude is 
overdamped by the critical layer compared to the classical 
case where only radiative and viscous dampings are taken 
into account; 

- in the second case, the fluid is unstable and turbulent and the 
critical layer acts as a secondary excitation region. Indeed, 
through over-reflection/transmission (when \R\ > 1 and 
\T\ > 1) energy is taken from the unstable shear that in- 
crease the amplitude of an incident IGW. Moreover, even 
in the case of simple reflection and transmission where 
\R\ < 1 and \T\ < 1, this demonstrate the existence of 
IGWs "tunneling" through unstable regions as identified by 
Sutherland & Yewchuk (2004) in laboratory experiments and 
by Brown & Sutherland (2007); Nault & Sutherland (2007) 
in Geophysics. 



Therefore, these mechanisms opens a new field of investigations 
concerning angular momentum tranport processes by IGWs in 
stellar interiors. 

Indeed, even if according to our first evolutionary calculations 
with STAREVOL, only the first stable regime exists in solar- 
type stars, we expect to find stars where the unstable regime 
and possible tunneling or over-reflection/transmission take 
place. Moreover, while the formalism presented in this work 
is general, several uncertainties remain. First, concerning the 
stochastic excitation of IGWs by convection, the model used in 
our evolutionnary code certainly underestimates the wave flux 
since it considers only the volumetric excitation in the bulk of 
the convective envelope while convective penetration should be 
also taken into account. This will influence the measured action 
of critical layers since it is proportional to the initial IGWs 
amplitude. Then, only retrograde waves are simulated here, 
considering that prograde ones are immediately damped and do 
not penetrate deeply in the radiation zone. This should normaly 
not afl'ect our results because the critical layers we detected are 
located in the deep radiation zone, but formal equations take 
both types of waves into account. 

The last point to bear in mind is that no latitudinal dependence 
for the angular velocity is considered here. We explained the rea- 
son of this choice in the introduction but one must not forget this 
approximation. Finally, since our goal is to get a complete and 
coherent picture of the transport of angular momentum in stellar 
radiation zones for every stellar type or evolutionary stage, it will 
be important to extend this work to the cases of gravito-inertial 
waves, where the action of Coriolis and centrifugal accelerations 
is considered (e.g. Lee & Saio 1997; Dintrans & Rieutord 2000; 
iviathis 2009; Ballot et al. 2010) and of magneto-gravito-inertial 
waves, when radiation zones are magnetized (e.g. Rudraiah & 
Venkatachalappa 1972; Kim & MacGregor 2003; MacGregor & 
Rogers 201 1; Mathis & de Brye 2012). 
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Appendix A: Validity of the JWKB approximation 

The form of the equation to solve is: 



Appendix B: Mathematical treatment for the non- 
perfect fluid case 

Let us introduce 

Xlm = r e'^'v{t)du (102) 



where a and b are the limits of a domain which will be defined 
later. The equation of propagation can then be written 



1 



and 



0. 



(104) 



— =/(r)^(r). 



(93) 



The first step is to introduce the Liouville transformation (e.g. 
Olver 1974), i.e.: 



iy(r) = //4vF, and^(r) 



We deduce 

dW _ 1 
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and Eq. (93) becomes 
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In the present case we have /(r) = -ky(r) = 
The WKBJ approximation is available 
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where Ri is the Richardson number of the fluid defined in Eq.(3). 

At last, we obtain that the condition for applying the WKBJ ap- 

1 /(/ -h 1) 

proximation is |0| ^ 1 which leads to Ri » — . 

4 



The new variable u defined by v = ^ 
inal equation into 
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(97) To obtain the final equation, we define: 
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and we get the following Whittaker equation: 
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